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ABSTRACT 

The far-wing line shape theory within the binary collision and quasistatic 
framework has been developed using the coordinate representation. Within this formalism, 
the main computational task is the evaluation of multidimensional integrals whose 
variables are the orientational angles needed to specify the initial and final positions of the 
system during transition processes. Using standard methods, one is able to evaluate the 
7-dimensional integrations required for linear molecular systems, or the 7-dimensional 
integrations for more complicated asymmetric— top (or symmetric— top) molecular systems 
whose interaction potential contains cyclic coordinates. In order to obviate this latter 
restriction on the form of the interaction potential, a Monte Carlo method is used to 
evaluate the 9— dimensional integrations required for systems consisting of one 
asymmetric-top (or symmetric-top) and one linear molecule, such as H 2 0 - N 2 . 

Combined with techniques developed previously to deal with sophisticated potential 
models, one is able to implement realistic potentials for these systems and derive accurate, 
converged results for the far-wing line shapes and the corresponding absorption 
coefficients. Conversely, comparison of the far-wing absorption with experimental data 
can serve as a sensitive diagnostic tool in order to obtain detailed information on the 
short-range anisotropic dependence of interaction potentials. 
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I. INTRODUCTION 

In two recent papers, 1 ’ 2 we reported progress made to overcome the two main 
uncertainties remaining in the calculation of accurate far— wing line shapes. In the first 
paper, 1 we introduced the idea of formulating the theory using the coordinate 
representation instead of the state representation. The motivation resulted from efforts to 
calculate converged line shapes for complicated molecular systems. In previous 
calculations,*’ 4 the states of the two interacting molecules were chosen as the complete set 
of basis functions in Hilbert space. These basis functions are characterized by well defined 
values of the energy, but not locations in space. As a result, one had to diagonalize 
anisotropic potential matrices because the latter are off-diagonal in this representation. 
However, it is well known that to diagonalize a matrix exhausts computer resources very 
quickly as the size of the matrix increases. Despite attempts to arrange the matrix in 
block— diagonal form which reduces the calculational cost significantly, the goal of 
obtaining converged results within reasonable computer CPU and memory constraints 
could only be achieved for the simplest systems such as C0 2 — Ar. By contrast, using the 
coordinate representation whose basis functions are characterized by the definite locations 
in space but not values of the energy, the diagonalization of the anisotropic potential 
matrix becomes unnecessary and the main computational task is transformed to the 
carrying out of multidimensional integrations over the continuous variables needed to 
specify the initial and final orientations of the two interacting molecules. For systems 
consisting of two linear molecules, or one linear and one more complicated 
(asymmetric— top or symmetric— top) molecule, or two complicated ones, the dimensionality 
of the integrations is seven, nine, and eleven, respectively. Then, as many states as desired 
can be included in the calculations and the convergence criterion is transformed to the 
feasibility of calculating these integrations within reasonable CPU time. In comparison 
with the previous requirement,*’ 4 this new criterion is much less stringent. We applied this 
formalism for linear molecular systems (C0 2 — C0 2 and C0 2 — N 2 ) and showed that one 
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could evaluate the 7— dimensional integrations using standard methods and obtained 
converged results with modest CPU time. 1 

The other uncertainty is related to forms of the interaction potential that could 
incorporated in the theory. It is well known that the theoretical line shape depends 
sensitively on the potential, especially on its short-range part. However, in previous 
formulations because of the necessity of finding roots of the radial part, the potential 
models that could be treated were limited to those having simple forms and one would not 
expect them to model the detailed interactions well. Therefore, an extension to more 
sophisticated potentials was necessary. In order to achieve this, we developed an 
interpolation procedure in which one does not have to find the roots directly; this advance 
enabled us to implement more realistic interaction potentials. The details of this technique 
were given in the second paper, 2 and combining this improvement with the coordinate 
representation, allowed us to calculate converged results for far— wing line shapes of linear 
molecular systems using realistic interaction potentials. More specifically, in a related 
paper, 5 we showed that using a sophisticated potential, one could calculate absorption 
coefficients that were in good agreement with experimental data over a range of 
temperatures for C0 2 — C0 2 . 5 

However, if we extend our method to other complicated cases (e.g., H 2 0 — N 2 or 
HjO — H 2 0), we found that the calculation of the higher-dimensional integrations is not 
feasible using standard methods. This implies that for these systems using standard 
methods, one has somehow to reduce the dimensionality of integrations required in 
numerical evaluations or to explore alternative methods to evaluate the 
higher— dimensional integrals. 

In the second paper, 2 we presented progress made according to the first procedure. 
We found that for cases where the interaction potential considered contains cyclic 
coordinates, the integration over them can be carried out analytically, with the 
concomitant reduction of the dimensionality of the remaining numerical integrations. The 
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later can then be carried out using standard methods. Because the leading long-range 
anisotropic interaction contains such cyclic coordinates, one can adopt a short-range 
interaction having the same symmetry. In comparison with results computed with less 
sophisticated potentials for the same systems, the latter results are in better agreement 
with experiment. 

In the present paper, we report results made by an alternative calculation of the 
multidimensional integrations using Monte Carlo methods. As a first step, we focus on 
systems containing one linear molecule. We show that by using the Monte Carlo 
technique, one can calculate the required 9-dimensional integrations and obtain converged 
results for the far— wing line shape using computational resources comparable to those 
needed to compute the 7— dimensional integrations using standard methods. Given the fact 
that we can consider potentials without cyclic coordinates, the advantages of the present 
method are obvious. 

The organization of the present paper is as follows: In Sect. II we give a summary of 
the general formalism. Starting in Sect. II A, we outline the relation between the spectral 
density, the Fourier transform of the dipole moment autocorrelation function, and the 
band— averaged line shape. In Sect. II B, we discuss the density matrix expressed in the 
coordinate representation that is central to the calculation of the far— wing line shape; we 
derive explicit expressions applicable for an asymmetric— top molecule and calculate the 
density matrix of H 2 0 for several different temperatures. We present representative 
two-dimensional plots in order to illustrate the main contour features. In Sect. II C, we 
discuss in detail the interaction potential used in the present calculation without the 
restriction of cyclic coordinates. This consists of the long-range dipole— quadrupole 
interaction varying as r' 4 , where r is the separation between centers of mass of the 
interacting molecules, an attractive dispersion term varying as r' 0 containing an adjustable 
strength, and a short-range exponential site— site repulsive interaction containing four 
adjustable parameters. In Sect. II D, we discuss the Monte Carlo calculation of the 
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9— dimensional integrations and the results obtained. In the final Sect., we present a brief 
discussion of these results and the conclusions drawn from the present study, including 
reasons why previous interaction potentials utilized for the calculation of the Lorentzian 
half— widths 8 ’ 7 are not suitable for the far— wing calculation. In fact, the accurate 
representation of the far— wing line shape places severe restrictions on the short-range 
interaction and can be used to discriminate between various potential models giving 
comparable results for other calculated quantities. 

H. THE GENERAL FORMALISM 

A. The spectral density and the band— average line shape function 

As shown in a previous study, 4 based on the binary collision approximation valid in 
the low— density limit, one is able to derive an expression for the spectral density F(u;) as a 
summation of contributions from individual lines with the common line shape (i.e., the 
band— average line shape) denoted by X(u): 

F <“> = ? I e — *-»— M ‘i)/ 2kI { (a , i ), /(*-<■.„) 

o;ij >o J 

+ (u> + W ..)5' , '( £J+ (!) 
where the indices i or j represent all the quantum numbers necessary to specify the energy 
levels of the absorber molecule, are the resonance frequencies, p i are the diagonal 
elements of the density matrix, and fi » are the reduced dipole matrix elements. Within the 
quasistatic approximation valid for the far— wing line shape, one can calculate the 
band— average line shape X(u) from knowledge of the interaction between two molecules. 

By choosing the Z axis of the space— fixed frame along the separation between the two 
molecules, the interaction V(r, fl a , fl b ) depends on the orientations of the molecules 
represented by fl a and fi b , respectively, and on r. Within the quasistatic approximation, 
the translational motion is treated classically while the internal degrees of freedom are 
treated quantum mechanically. In Hilbert space associated with the internal degrees, one 
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denotes the eigenvalues and eigenvectors of V(r, ft a , fl b ) by and | (>, respectively; thus 

V(r,n„n b )K>=G t (r)K>, (2) 

where r appears in the eigenvalues G^(r) as a parameter. Then, one is able to derive an 
explicit expression for the band— average line shape X(J): 

}{u) = 47T 2 n b or* \<(Wp^n\f}>\ 2 


* r 


2 

c 




(v 


q — [ G <,( r c)+ G t)( r c)]/2kT — Viso(rc)/kT y ^ p | 


ij 



(3) 


where n b is the number density of the bath molecules, G^(r) denotes gjlG;( r ) — G^r)], 
and r c are roots of the equation 


G;(r c ) - G^rJ = u. (4) 

In order to carry out numerical calculations of X(u) from Eq. (4), one first has to 
find the eigenvectors and eigenvalues of the potential and then sum the contributions 
resulting from all possible combinations of these eigenvectors representing the initial and 
final positions of the system during the transition processes. It has been shown that if the 
states of the system are chosen as the complete set of basis functions in Hilbert space, one 
is unable to obtain converged results except for the simplest systems such as C0 2 — Ar 
because a diagonalization procedure required exhausts computer resources quickly. In 
order to overcome this difficulty, a new method was developed recently in which the 
eigenfunctions of the orientations of the system are chosen as the complete set in Hilbert 
space. 1 From the theoretical point of view, these two methods are equivalent since they 
differ only by the choice of representation, but in practical calculations, the latter has a big 
advantage over the former. In fact, with the new method the diagonalization procedure is 
unnecessary and the main computational task is transformed to the carrying out of 
multidimensional integrations over the continuous orientational variables. With this 
formulation, we have shown that using standard methods, one is able to calculate 
converged results for systems consisting of two linear molecules in which the dimensionality 
of the integrations is seven. 1 ’ 5 For systems involving more complicated molecules, such as 
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H a 0 — Nj and H 2 0 — H 2 0, the dimensionality of the integrations becomes nine and eleven, 
respectively, since more angular variables are needed to specify their orientations. Usually, 
standard methods are not applicable in order to evaluate such high multidimensional 
integrations. Recently, we have shown that in cases where the interaction potentials 
contain cyclic coordinates, one can carry out the integrations of the density matrix over 
these coordinates analytically first and obtain the "averaged" density matrix. 2 Then, with 
the "averaged" density matrix, one is able to reduce the dimensionality of the remaining 
integrations for both H 2 0 — N 2 and H 2 0 — H 2 0 to seven, the same as those for linear 
molecular systems. As a result, one can obtain converged results using the standard 
methods. Although the leading terms of the long-range interactions of H 2 0 — N 2 and H 2 0 — 
H 2 0 do have such cyclic coordinates, to require the other parts, especially the short-range 
part which is not known well at present, to have the same symmetry is a simplifying 
assumption which one would like to avoid. 

It has been known for a long time that Monte Carlo methods can be applied for 
evaluating high multidimensional integrations. Given the fact that a calculated line shape 
ranges from line center to more than 1000 cm* 1 away and requires a lot of points to depict 
well, we note that the complexity associated with the line shape calculation is greater than 
the evaluation of a few 9— or 11— dimensional integrations. As a first attempt, without 
introducing any assumptions on the symmetry of the H 2 0 — N 2 interaction, we use the 
Monte Carlo method to directly evaluate the 9-dimensional integrations and calculate the 
line shape. By doing so, instead of the "averaged" density matrix of H 2 0, the complete 
density matrix has to be used. 

B. The density matrix in the coordinate representation 

Within the band— average approximation, we have shown that one can replace 
\<<\4PJ>\A 1 2 by |<CIV?a?bl V> 1 2 i n the expression for the band— average line shape. 1 

In later derivations, formulas without the subscripts a or b attached are understood to be 
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applicable for both the absorber and bath molecules. 

In the following derivations, we adopt the simplifying notation that is used to 
represent the orientation of the molecule of interest. For N 2 , corresponds to 0^ and ip^. 
Alternatively, for the asymmetric— top H 2 0 molecule, it corresponds to a^, /?<-, and 7 ^. 
Then, expressing | (> as 


IO = I «(«, - n .0 «( Q b - n H )>. (5) 

one is able to write I <CIV/’ a t'bl b> 1 2 as the product of two factors 

l<CIVvbh>l , = l<£(n a - n. t ) | JFJ «(n a - «,,)> | ■ » 


« I < *(«b - Ob() I «(»b - n bs)> I ’• (6) 

The explicit expression for each factor depends on which class the molecule of interest 

belongs to. For linear molecules, we have shown that 1 

I «(n - n ; )| v?| i(i> - n„)> | 1 = l a l pjcose, w ), (7) 

L 


where P L (cos0 ( ^ ) are Legendre polynomials with L = 0 , 1 , 2 , • • •, and 0^, is the angle 
between and fl^, viz. 

cos©^, = cos^ cos^ + sinfy- sinfl^ cos(y^ — ( 8 ) 

The coefficients A are given by 

L 


A l = TO 57 J l (2j + l)(2j> + 1) ^ e (E(j)+E(j’)l/2kT c ’ (j j> Lj 0 0 0 ), 


(9) 


JJ 

where Q is the partition function, g is the nuclear spin degeneracy factor, E(j) are the 
energies of the state corresponding to the angular quantum number j, and C(j j’ L, 0 0 0 ) is 
a Clebsch-Gordan coefficient. The absolute square of the density matrix for N 2 calculated 
from Eqs. (7) and (9) at T = 296, 338, and 430 K with a cut-off = 80 are shown in 
Fig. 1 . As shown in the figure, there are two sharp and uneven peaks located around 0 ( a , 
= 0 and 180 degrees, respectively, and the magnitude decreases very fast as 0 ( differs 
from these values. By comparing the profiles obtained from different temperatures, it is 
obvious that the higher the temperature, the higher and narrower the peaks become. The 
more detailed discussions were reported previously . 1 
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For a symmetric— top molecule whose wavefunctions are given by 8 

| jkm> = (- l) m - k Di m _ k (a, 0 , 7 ), (10) 

where k (o, 0 , 7) (= e “ k (j 9 )e"^ k ^) is the rotational matrix, we have shown 2 

that 

— — n ir) )>| 2 = ^ A LK 7(^7,)). (11) 

LK 

In Eq. (11), a ( < - T)) , f 3 { ( - T)) , and 7 ( < - T)) are the three Euler angles used to represent a rotation 
resulting from two successive rotations, i.e., 

®-( a ( Ctj) > ;t)) > ;ti> ) = &(,’ 1 <) ^(°fy ^t|> 7t))> (12) 

the summation index L = 0, 1, 2, • • • and the index K = — L, — L + 1, • • • , L — 1 , L, 

respectively, and the coefficients A are given by 

LK 

A LK = MFQ II + V ^hh-K 

* e - (E(j,k)+E(j , ,k-K)]/ 2 k b T c 2 q L> k K _ k (13) 

where E(j,k) are the energies of the state labeled by the quantum numbers j and k, and g k 
is its nuclear spin degeneracy factor. We note that when the quantum number k and the 
Boltzman constant k appear in same expressions, a subscript b is attached to the latter. 


Similar expressions can also be derived for the asymmetric— top molecule whose 

wavefunctions are given in terms of an expansion of symmetric— top wavefunctions, 9 

| jrm> = J U kT | jkm>. ( 14 ) 

k 

In this case, we have shown 2 that 


| <«(n - n t ) |y£| 6(0 - a n )> I ! = n A M- D i,i'(“cw^wT, w ). (15) 

L KK’ 

where ^ 7 ( ;ti> **6 defined by Eq. (12) and the summation index L = 0, 1, 

2 , • • • , and both indices K and K’ run from — L to L. In the above expression, the 
coefficients A^ k , are given by 


L _ 1 
KK’ ~ SiFQ - 


l l (-l) K + K ’(2j + l)(2j' + l)^e- 


[E(j,r)+E(j’,r’)]/ 2 k b T 


. . > 

JJ 


rr’ 
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I(-l> 


k + k’u J TT ^ TT J ' 

U kr U k’r U k’-K’ 


r’ U k-K r 


, C(j j’ L, k K-k K) C(j j’ L, k’ K’-k’ K’), (16) 


kk’ 

where E(j,r) are the energies of the state, and is its nuclear spin degeneracy factor. By 
comparing Eqs. (7), ( 11 ), and (15), i.e., the expressions for the factors associated with the 
linear, the symmetric— top and asymmetric— top molecules, respectively, it is obvious that 
the more complex the molecule, the more complex the expression becomes. In addition, for 
the linear molecule, the absolute square of the density matrix represents a one-dimensional 
distribution over the angle © ( . Meanwhile, for the symmetric— top and asymmetric— top 

molecules, it can be understood as a three-dimensional distribution over the Euler angles 
a ( , P { ^ , and 7 { . For later convenience, we use the simple notation F(a, P, 7 ) to 

represent it. We note that the normalization requirement is given by 


r27T r7T r 2x , 

Jo J 0 j 0 F ( q > A 7 ) sinP da dP dy = (17) 

It has been shown 1 ’ 2 that one can easily include as many states as desired in 
numerical calculations for the density matrix of the linear molecules C0 2 and N 2 and also 
for the "averaged" density matrix of H 2 0. In the present study, we note that one can still 
include all the states available in the calculations with affordable CPU costs. Although in 
order to calculate the density matrix for complex molecules there are many more individual 
coefficients A and A^ , needed, and their expressions become more complicated as shown 
by Eqs. (13) and (16), respectively, the calculation of these coefficients is straightforward. 
In addition, one does not need to calculate all of them since some of them are identical and 
others are zero. For the H 2 0 molecule of interest in the present study, it turns out that 
due to the symmetry of u£ r all the coefficients A^ k , are zero unless their indices K and K’ 
have the same evenness or oddness. In addition, for the non-zero coefficients the following 


symmetries are valid: 


A l — A l — A l 
KK’ K’K — K-K’ ■ 


(18) 


As a result, if j max is the highest angular quantum number of the states included in the 
calculation, one only needs to calculate the coefficients A^ k , for L = 0, 1 , 2, • • • where 
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L max = 2 jmax> for K = 0, 1, 2, • • • L, and for K’ = - K, - K + 2, • • • , K - 2, K. The 
number of distinct coefficients is equal to (L max + 1)[1 + L max (L max + 5)/6]. As an 
example, for the pure rotational band of H 2 0, the highest angular quantum number of the 
initial states listed in the HITRAN 92 database 10 is 23. If one includes all these states, 
there are 18424 coefficients needed to be evaluated. In comparison with the calculation of 
the "averaged" density matrix in which there are only 47 (= 2j max + 1) coefficients A L 
needed, of course, this requires a lot of more computation requiring several dozen hours to 
complete using a workstation. But, fortunately, for a specified temperature it needs to be 
done only once. 

When all these A^ k , are available, one can easily obtain the density matrix given in 
the coordinative representation since the the expression for the rotational matrix is known. 
The explicit expression for F(a, 0, 7) used in the numerical calculations is given by 
F (“. A 7) = n l hr C0S K K + K ’)(“ + 7)/2] cos[(K - K')(q - 7)/2] (19) 

LH’ 

where e KK , = 1 for K = K’ = 0; e KK , = 2 for |K’| = K; Cj[K , = 4 for |K’| * K, and the 
ranges of the indices L, K, and K’ are from 0 to L,^, from 0 to L, and from — K, — K + 2, 

• • • , K — 2, K, respectively. We note that in this expression, the Euler angles a and 7 
appear not by themselves, but in the combinations (a + y)/2 and (a — i)/2. In addition, 
(a + 7)/2 is associated with K + K’, and (a - 7) /2 is associated with K — K’. More 
discussion about this point is given later. 

It has been shown 2 that the profile of the "averaged" density matrix of H a O which is 
a one-dimensional distribution to some extent looks like a 6 function. One expects that 
the profile of the density matrix of H 2 0 itself which is a three-dimensional distribution 
exhibits more detailed features. Of course, it is impossible to plot F(a, 0, 7) in one figure 
and one has to present it using a series of two-dimensional contours with one argument 
fixed. 

With Eq. (19), we calculate several two-dimensional distributions F(a, 0 O , 7) over 
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the Euler angles a and 7 at T = 296 K obtained with 0 O = 0, 20, 40, and 60 degrees, 

respectively, and present their corresponding contour plots in Fig. 2. From this figure, one 

can draw two conclusions. By comparing these plots, one can easily find that the 

magnitudes of the density matrix decrease very fast as 0 O increases. This is fully expected 

since the profile of the "averaged'' density matrix known from a previous study exhibits 

this behavior. In fact, by carrying out a two-dimensional integration of the density matrix 

over the Euler angles a and 7, one expects to obtain a one-dimensional distribution of the 

density matrix over the remaining angle 0. This results in exactly the "averaged" density 

matrix as easily shown by 
p2tt r 2 w 

j Q J Q F(a, 0 , 7) da d 7 = 4^ J A L QQ Dj 0 (0, 0, 0) 

L 

= l A l P l (cos/?). (20) 

L 

In deriving the last step of Eq. (20), the fact that 4 x 2 Aqq = A L has been used and the 
explicit expression for A was given in our previous study. 2 Secondly, by analyzing each of 
these contours, it is obvious that these F(a, 0 O , 7) exhibit symmetry with respect to the 
axes (a + 7)/2 and (a — 7) /2. It is interesting to note that if one introduces the two 
arguments u = (a + 7)/2 and v = (a — 7)/2, the magnitudes of F(a, 0 O , 7) depend very 
sensitively on the former, but are insensitive to the latter, especially for the cases where 0 O 
is small. This implies that instead of a and 7, it is better to choose u and v as the two 
arguments. The same conclusion is also suggested by Eq. (19) as mentioned above. We 
use the symbol G(0, u, v) to represent the distribution of the density matrix over 0, u, and 
v. We note that the values of u and v vary in a lozenge-shaped area, the former’s 
maximum and minimum equal to 0 and 2x, respectively, and the latter’s are — x and x, 
respectively. 

Since G(/3, u, v) has one insensitive argument v, its profile can be well represented 
by one contour plot of TS(0, u) resulting from averaging over v. Due to the symmetry TJ(0, 
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2x — u) = Tj(0, u), one can calculate values of u) only in the range of v from 0 to 7 r. 

In this case, u) is defined by 

u) = 25 f G(/3, u, v) dv. (21) 

With Eq. (19), one is able to carry out the above integration analytically and derive an 
explicit expression for G(/?, u) but this is not presented here. We calculated G(/7, u) for T 
= 296 K, 338 K, and 430 K and present the results obtained in Figs. 3-5, respectively. In 
order to visualize the complete profiles of the density matrix, one can imagine that except 
for small deviations which can be ignored in small /? areas, the profiles shown in Figs. 3 — 5 
extend along another dimension, i.e., the v axis which is perpendicular to the /? — u plane 
and is missing in these figures. Furthermore, the range of the extension along the v axis is 
equal to [— 1 u, u], which varies from the minimum 0 at u = 0 to the maximum 2x at u = x. 
From these figures, it is obvious that the density matrix of H 2 0 exhibits more structures 
than that shown by the "averaged" density matrix. Along the u axis, there are 5 sharp 
peaks located at (0, 0), (0, x/2), (0, x), (0, 3x/2) and (0, 2x), respectively. The largest is 
the one at (0, t). This conclusion results not only from its apparent size in the plots, but 
also from the fact that it extends along the v axis more than the others. With respect to 
its distributions along the /? axis, in general, the density matrix of H 2 0 decreases very fast 
as P increases which has been shown previously for the "averaged" density matrix. On the 
other hand, by comparing the results obtained from different temperatures, one concludes 
that as the temperature increases, the peaks become sharper, i.e., their heights are higher 
and their widths are narrower. We note that in the contour plot obtained at T = 430 K, 
the lowest scale is reduced from 1.0 * 10‘ 1S used for T = 296 K and 338 K to 1.0 * 10' 10 . 
Otherwise, calculation errors add some meaningless structures in areas where the 
magnitudes are below this limit; this indicates a general trend that for higher temperatures 
more states with higher j values have to be added in the calculations. 

In order to carry out line shape calculations, one can calculate G(/?, u, v) and store 
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the results in files. The areas with high /? values can be excluded from consideration. One 
can further reduce the ranges of u and v due to symmetry considerations. One has to 
choose high resolutions for and u, but a lower one for v. The details are not presented 
here for brevity. 

C. More general potential models 

In our previous study 2 , we have developed techniques to deal with complicated 
interaction potential models, specifically those for which their dependence on r and on fi 

& 

and Hjj interweave, and the dependence on r takes a complicated form. For example, we 
can adopt a flexible site-site model, 

V sr (r, n a , n b ) = l l V„(i„), (22) 

iea jeb 

in modeling the short-range part of the anisotropic interaction. In this site— site model, 
the indices i and j run over force centers of the absorber molecule a and the perturber b, 
respectively, r^ is the distance between the center i and the center j, and V ij (r i j) represent 
individual interactions between them. The flexibility of the site-site model can be 
achieved by choosing different numbers of the centers, different locations of the centers, 
and different forms for Vy. With respect to the long-range part of the anisotropic 
interaction, one usually adopts the leading term of the multipole expansion which, in the 
present case, is the dipole-quadrupole interaction given by 

v dq( r > n a> fi b) = 5^? [ C0S Ai ( 3 cos 2 0 b - 1) - 2 sin/? a sinfl b cos0 b cos(a a - y^)], (23) 
where is the dipole moment of H 2 0 lying along the Z axis of the body-fixed frame, and Q 
is the quadrupole moment of N 2 . 

It has been shown that in order to make multidimensional integrations tractable 
with standard methods, we had to limit potential models to those containing cyclic 
coordinates. The dipole-quadrupole interaction of H 2 0 — N 2 has a cyclic coordinate 7 
By assuming that in the site-site model the force centers of H 2 0 are located along the Z 
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axis of the body— fixed frame, V sr (r, fl a , fl b ) has the same symmetry as V dq (r, fl a , fl b ). 
With this restriction, one can introduce the "averaged" density matrix and reduce the 
original 9— dimensional integrations to 7— dimensional ones. Then, standard methods 
become applicable to carry out the integrations. 

Since at present one does not know the short-range anisotropic interaction of H 2 0 — 
N 2 well, it is desirable to consider more general models. Therefore, we adopt a site-site 
model having three force centers for H 2 0: two located at the two H atoms and one at the 0 
atom. With respect to the N 3 molecule, we assume that there are two force centers and 
they are located on the two N atoms, the same as we did previously. Concerning the form 
of V i j(r i j), we assume that V i j(r i j) = Ajj e — to represent the repulsive interaction 
between force centers where Ajj and p i j are adjustable parameters. In order to take into 
account attractive dispersion interaction, we include an isotropic term given by — B/r 6 
where B is an adjustable parameter. As a result, the potential V(r, fl a , fl b ) considered is 
given by 

V(r, <1 S , n b ) = l l Ay e- + v dq (r, n„ n b ) - f, , (24) 

iea jeb 

and it contains 5 adjustable parameters: A UH , p uv , A n „, p nv and B. The dipole moment of 
H 2 0 and the quadrupole moment of N 2 are well known and the values used in the present 
calculations are 1.8546 D and 1.40 D A, respectively. Meanwhile, the bond angle of H 2 0 is 
104.52 degrees and the separation between the 0 atom and the H atom is 0.9572 A; 11 the 
separation of two N atoms of the N 2 molecule is 1.10 A. 12 By making these choices, V(r, 

O a , ft b ) depends on all 5 angular variables used to specify the orientations of H 2 0 — N 2 . 

This implies that V(r, ft a , fl b ) does not contain any cyclic coordinates and it is necessary 
to evaluate the 9— dimensional integration. 

D. A Monte Carlo calculation of the line shapes 

By choosing the eigenfunctions of the orientations of H 2 0 — N 2 as the complete set 
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of basis functions, the summation over ( and rj in Eq. (3) becomes 9— dimensional 
integrations over 0 a; , 7 a(; , 0 b <-, a at) , 0 aT) , 7 aT ,, 0 bT ,, and ^ in which the first four 
(including a a( - = 0) specify the initial orientations of H 2 0 — N 2 and the last five specify the 
final ones. We note that due to the rotational symmetry of the whole system, one can 
always assume o a ^ = 0. 

It has been known for years that the Monte Carlo method is an alternative way to 
calculate multidimensional integrations, especially in cases where their dimensionalities are 
so high that standards methods do not yield reliable results. We rely on an algorithm 
called VEGAS which is widely used in elementary particle physics. 13 As a Monte Carlo 
algorithm, VEGAS is primarily based on important sampling, but also does some stratified 
sampling. The important sampling is crucial in the present calculations since the 
distributions of the integrand over its 9 integration variables are extremely nonuniform. In 
fact, as shown by Eq. (3) it contains the factor | <C| V , P a / J b/ 1 l *?> 1 2 which is the product of 
the density matrix of H 2 0 and the density matrix of N 2 . As shown by Figs. 1 — 5, these 
density matrices are highly concentrated in narrow areas. As a result, one must use the 
important sampling to achieve variance reduction in Monte Carlo computations. This can 
easily be done by setting input flags in VEGAS. With VEGAS, the subroutine vegas 
performs m (the flag 1) statistically independent evaluations of the desired integral, each 
with N (the flag 2) integrand evaluations. While statistically independent, these iterations 
do assist each other, since each one is used to refine the sampling grid for the next one. 

The results of all integrations are combined into a single best answer and its estimated 
error. Interested readers can find its basic features in the book of Numerical Recipes in 
Fortran 77. 13 

We have made some modifications necessary in VEGAS resulting from special 
requirements in the line shape calculation. As mentioned previously, the complexity 
associated with the line shape calculations is far beyond the need to evaluate a few values 
of the 9— dimensional integrations In practice, since the line shape is a function of 
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frequency, one has to use a lot of points to depict it well. Therefore, the number of the 
9— dimensional integrations to be evaluated simultaneously by the algorithm is around 100 
or so. 

In order to verify the applicability of Monte Carlo algorithm in the line shape 
calculation, we have made several tests in which the potentials considered contain cyclic 
coordinates and 9-dimensional integrations can be reduced to 7— dimensional ones. In 
these cases, we use two methods to calculate the line shapes. First, by introducing the 
"averaged" density matrix and using standard methods, we carried out the reduced 
7— dimensional integrations and obtained converged results. On the other hand, we used 
the Monte Carlo algorithm directly to evaluate 9-dimensional integrations. It turns out 
that in general, by taking around 10 6 or 10 7 random selections of a set of the 9 
orientational variables, the Monte Carlo algorithm can produce line shapes without 
significant differences from those obtained before. 2 The CPU times required are modest on 
a workstation. 

Then we applied the Monte Carlo method to cases where the potentials do not 
contain any cyclic coordinates. In order to verify the applicability, one can compare results 
obtained from the Monte Carlo algorithm with different random selection numbers and 
check whether results have converged. We present in Fig. 6 two line shapes obtained at T 
= 296 K with the same potential parameters but with different numbers of the random 
selections. The first one is obtained from 3.2 * 10® random evaluations of the integrand, 
i.e., a call of VEGAS with m = 3 and N = 4 » 10 5 immediately followed by another call 
with m = 1 and N = 2 * 10®. The second is from 9.2 « 10® random evaluations using a call 
of VEGAS with m = 3 and N = 4 * 10® immediately followed by another call with m = 1 
and N = 8 * 10®. From Fig. 6, it is obvious that there are only negligible differences 
between these two line shapes. Therefore, one can conclude that using approximately 10 7 
random selections of the set of variables, the calculated line shapes are converged. It is 
interesting to note that the same number of terms is required to evaluate 7— dimensional 
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integrations using standard methods. This implies that with comparable computer 
resources (several hours CPU time at a workstation) required previously for 7— dimensional 
integrations, the Monte Carl algorithm enables us to evaluate 9-dimensional ones. 

After making convergence tests, we know how many random evaluations are needed 
in order to obtain reliable results. Thus the only difficulty left in the line shape 
calculations, and consequently in the absorption coefficient calculations, is the choice of a 
realistic potential model. We adopt the potential form given by (24) and begin to search a 
set of potential parameters that enables us to get a good fit to the experimental results of 
Burch et al. 14 at T = 296 K. We find that by adopting A /k = 1.2 * 10® K, p = 0.28 
A, A QK /k = 1.5 « 10 8 K, p 0(| = 0.26 A, and B/k = 1.38 * 10® K, one is able to obtain 
satisfactory results. In addition, we note that because the short-range model has an 
isotropic contribution, which when combined with the isotropic dispersion interaction 
— B/r®, gives "effective" Lennard— Jones parameters of tr eff = 3.57 A, and e eff /k = 207 K, in 
reasonable agreement with other estimates. Based on this potential model, we calculated 
the the band— averaged line shapes of H 2 0 — N 2 for several temperatures T = 296 K, 338 K, 
and 430 K. The line shapes at T = 296 K plotted in Fig. 6 are from these calculated 
results. The line shapes at T = 338 K and 430 K, together with that at T = 296 K, 
obtained from a call of VEGAS with m = 3 and N = 4 * 10® immediately followed by 
another call with m = 1 and N = 4» 10®, are presented in Fig. 7. The N 2 — broadened 
absorption coefficients at T = 296 K in the spectral region 300 — 1100 cm' 1 based on 
HITRAN 92 data are plotted in Fig. 8, together with the experimental results of Burch et 
al. 14 Similarly, the calculated absorption at T = 430 K and comparison with experimental 
data 14 are presented in Fig. 9. As shown by these figures, the agreement with the 
laboratory values is reasonably good. 

III. DISCUSSIONS AND CONCLUSIONS 

In summary, there are three main theoretical improvements made recently in line 
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shape calculations. The first is that instead of the representation constructed from the 
states of the system, one uses the coordinate representation enabling one to avoid 
convergence errors resulting from the truncation of states included in calculations, and the 
failure to maintain the commutativity of the dipole moment and the anisotropic potential. 
We note that in practical calculations, the ability to overcome the convergence problem is 
crucial for the accurate calculation of the corresponding far— wing line shape. 

The second improvement is the interpolation procedure introduced which obviates 
the need to find the roots in Eq. (4) and enables one to handle complicated potential forms. 
Given the fact that the calculated far— wing line shapes depend sensitively on the potential, 
especially on its short-range anisotropic part, and that the latter can not be represented 
well by simple models consisting of the Lennard— Jones isotropic potential and the leading 
term of the multipole expansion, the ability to incorporate more sophisticated potential 
models is crucial for realistic calculations. For linear molecular systems, it has been shown 
that the main computational task is to carry out 7-dimensional integrations and the later 
can be evaluated using standard methods. However, for systems involving complicated 
molecules, the generality of potential models is limited to those having cyclic coordinates. 
This restriction permits one to introduce the "averaged" density matrix in order to reduce 
the dimensionality of integrations to be evaluated. Otherwise, standard methods are not 
applicable for carrying out the required 9— or 11— dimensional integrations. 

Finally, the Monte Carlo algorithm has been tested for line shape calculations. 

Based on numerical tests, one concludes that this method can be used, at least to evaluate 
9— dimensional integrations. As a result, one can utilize sophisticated potential models for 
systems consisting of one asymmetric— top (or symmetric— top) molecule and one linear 
molecule and obtain converged results for the far— wing line shapes and corresponding 
absorption coefficients. 

It has been shown in our previous study 2 that the far— wing shape is very sensitive to 
the angular gradients of the potential and, more specifically, to the detailed energy contour 
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area at which large angular gradients of the potential are exhibited while the potential 
values are relatively small or even negative. We expect that such features would not fully 
manifest their effects on other physical measurements (e.g., the second virial coefficient, 
the Lorentzian halfwidth and shift, etc.). This implies that the far— wing line shape can 
serve as a sensitive diagnostic tool to obtain detailed and desirable information about the 
interaction between two molecules. Given the fact that we can treat all potential models 
for H a O — N 2 and we have full confidence that the results obtained are the "true" values 
corresponding to these models, this role becomes even more valuable. 

In the present work, we have tested some other models, such as a site-site model 
whose expression is given by 

V(r, n„ n b )= l l 4c,j [(fl)“ - (fii)Y (25) 

iea jeb 1J 1J 

where the individual interactions between force centers take the Lennard— Jones form and 
Cy and (Ty are corresponding parameters. It turns out that with the same parameters used 
in Lorentzian halfwidth and shift calculations, 6 ’ 7 the magnitudes of the calculated far— wing 
line shape are too large and predict much more absorption when compared with the 
measurements. By analyzing contour profiles of this model in detail, we find that its 
angular gradients are very large at some areas where the potential values are relatively 
small. This results from the Lennard— Jones form. It is well known that around r = a, the 
value of the Lennard— Jones potential is close to zero, but its gradient in r is large. For 
each components Vy, the same feature appears around r^ = cr iy On the other hand, 
although the distance r between two molecules is assumed to be fixed during the transition 
process, all the distances between force center i and center j, i.e., r A j vary as the two 
molecules change their orientations. It happens that for some special positions of the 
system, one or more of these r^ are close to a iy Then, small rotations of the orientations 
could result in significant changes in the values of Vy. As a result, V(r, Q a , fi b ) could 
possess very large angular gradients and relatively small values near these special areas. 
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As mentioned above, most of significant contributions to the magnitudes of the 
far— wing line shape come from positions of the system such that around their energy 
contour area, large angular gradients of the potential are exhibited while the potential 
values are relatively small or even negative. The latter requirement can be easily 
understood since the system would have few statistical chances to occupy positions having 
large positive potential values. Therefore, it is not surprising that the potential model 
constructed from the Lennard— Jones form produces too large magnitudes of the far— wing 
line shape. The failure to predict the absorption indicates that the actual interaction is not 
well represented by this model. However, we note that this does not devalue its 
applicability for predicting different physical quantities which are insensitive to these 
features, but sensitive to other ones. The profile of the interaction potential between two 
molecules has a complicated nature which cannot be ascertained from only a single physical 
measurement; one has to combine all diagnostic means available to map its complexity. 

Although the work performed in the present paper is for H 2 0 — N 2 , similar studies 
can be carried out for other systems, such as H 2 0 — C0 2 . However, one challenge still 
remains, i.e., how to deal with the H 2 0 - H 2 0 system for which one needs to evaluate 
11— dimensional integrations if the potential considered does not contain any cyclic 
coordinates. 

In general, to extend the Monte Carlo computations to higher dimensions may be 
possible. It will require more random selections to achieve converged results, but not 100 
times more (a typical amount obtained by selecting 10 points along each additional 
dimension) as standard methods do. However, there is another difficulty faced by the 
Monte Caxlo computations resulting from the integrand itself, i.e., that its distributions 
become even more nonuniform. In fact, the factor | <ClV/ , a / , b/ i l 1 3 ^ or — 
distributes nonuniformly over 6 variables in comparison with the H 2 0 — N a case where it 
distributes nonuniformly only over 4 variables. 

In addition, it has been shown in our previous studies 1 ’ 2 that it is essential to 
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separate the sensitive variables of the integrand from insensitive ones since this separation 
enables one to choose different resolutions in evaluating multidimensional integrations 
using standard methods. The same conclusion is also true in the Monte Carlo 
computations. In this case, the separation enables one to tailor the important sampling 
and to reduce the variance dramatically. This separation is achieved by representing the 
final orientations of the system labeled by rj in terms of the body— fixed frames instead of 
the space— fixed frame. The body— fixed frames introduced here are those attached to the 
two molecules at their initial orientational positions. More explicitly, for the linear 
molecule, one chooses 6^, 6 ( ^ , and ^ T)) instead of 6 tp^ t 6^, and <p^, for the 

symmetric-top or asymmetric— top molecule, one chooses a^, /?<-, 7 ^, > , and 

7 ( ^, instead of a^, /?<-, 7 ^, a^, 0^, and 7 ^. However, it is necessary to carry out some 
algebra when one calculates potential values with these variables since the expressions for 
potentials are always given in terms of the space— fixed frame; thus, one has to rotate the 
body— fixed frame back to the space fixed one. The more complex the molecule, the more 
work is required. 

Based on the arguments given above, it is obvious that more random selections and 
more CPU time will be required to achieve converged line shapes for H 2 0 — H 2 0. 

However, without numerical tests, it is difficult to assert whether the Monte Carlo method 
is applicable for this system. We plan to pursue this problem in further research. 
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LIST OF FIGURES 

Fig. 1. The absolute square of the density matrix of N 2 given in the coordinate 

representation calculated at three temperatures for = 80: the solid line at T = 
296 K; the dashed - dotted line at T = 338 K; and the dotted line at T = 430 K. 

Fig. 2. The two-dimensional distributions of the density matrix of H 2 0 over the Euler 

angles a and 7 obtained at T = 296 K for j max = 23. The value of the Euler angle 0 
is fixed and there are four plots presented here corresponding to 0 = 0, 20, 40, and 
60 degrees, respectively. 

Fig. 3. The two-dimensional distribution of the density matrix of H 2 0 over the two 

sensitive angles 0 and u [= (0 + ?)/2] obtained at T = 296 K for j max = 23. This 
distribution results from averaging the density matrix of H 2 0 over the one 
insensitive angle v [= (a — The range of u is [0, 2r] and there is a symmetry 

respect to the axis u = t. The plot presented here covers 0 < u < x only. 

Fig. 4. The same as Figure 3, except T = 338 K. 

Fig. 5. The same as Figure 3, except T = 430 K. In comparison with that used in Figures 
3 and 4, the lowest scale is reduced to avoid meaningless structures which result 
from numerical errors. This indicates a general trend that for higher temperatures 
more states with higher J values have to be added in the calculations. 

Fig. 6. The N 2 -broadened far-wing line shape of H 2 0 (in units of cm' 1 atm" 1 ) as a 
function of frequency u> (in units of cm' 1 ) for T = 296 K calculated with two 
different numbers of random selections in the Monte Carlo calculation; they are 
obtained from 9.2 * 10 6 and 3.2 * 10® random evaluations of the integrand and are 
represented by the solid and dashed curves, respectively. 

Fir. 7. The N 2 -broadened far-wing line shape of H 2 0 (in units of cm' 1 atm' 1 ) as a 

function of frequency u (in units of cm' 1 ) for T = 296 K, 338 K, and 430 K; these 
are represented by the solid, dashed, and dotted lines, respectively. 

Fig. 8. The calculated N 2 — broadened absorption coefficient (in units of cm 2 molecule' 1 
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atm' 1 ) at T = 296 K in the 300 — 1100 cm' 1 spectral region is represented by a. For 
comparison, the experimental values are denoted by +. 

Fig. 9. The same as Fig. 8, except that T = 430 K. 
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